Phase 1: Exploration
Purpose: Systematic exploration of
senescence-associated gene lists from multiple sources, comparing them
with high-confidence DEGs from current study.
Gene list sources:
SAG orthologs from Zhang et al. 2014
Natural senescence genes from Zhang et al. 2014
SAG GWAS hits from Sekhon et al. 2019
Stay-green network from Sekhon et al. 2019
SAG hub genes from Ojeda et al. 2026
Photosynthesis genes from Ojeda et al. 2026
GO-based photosynthesis and leaf senescence gene sets, Fattel et
al 2024
# Load cross-reference tables for gene ID conversion
v3_v5 <- read.table(
file = file.path(paths$data, "B73v3_to_B73v5.tsv"),
sep = "\t",
header = FALSE
) %>%
rename(v3 = "V1", v5 = "V2") %>%
separate_longer_delim(cols = v5, delim = ",")
v4_v5 <- read.table(
file = file.path(paths$data, "B73v4_to_B73v5.tsv"),
sep = "\t",
header = FALSE
) %>%
rename(v4 = "V1", v5 = "V2") %>%
separate_longer_delim(cols = v5, delim = ",")
SAG Orthologs - Zhang et al. 2014
Purpose: Load SAG orthologs identified through
comparative genomics and convert to v5 gene IDs.
sag_orthologs_Zhang2014 <- read.csv(
file.path(paths$data, "SAG_orthologs.csv")
) %>%
janitor::clean_names() %>%
inner_join(v3_v5, by = c(gene_id_maize = "v3"))
cat("Total SAG orthologs (Zhang 2014):",
length(unique(sag_orthologs_Zhang2014$v5)), "\n")
## Total SAG orthologs (Zhang 2014): 47
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% sag_orthologs_Zhang2014$v5
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 15 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Upregulated Zm00001eb1… <NA> <NA> 2.28 4.18e- 8
## 2 -P Upregulated Zm00001eb3… <NA> Epimerase … 3.02 4.72e- 5
## 3 -P Upregulated Zm00001eb2… <NA> ATP-depend… 2.23 2.33e- 4
## 4 -P Upregulated Zm00001eb0… pco102923 Aspergillu… 3.68 2.14e- 3
## 5 -P Upregulated Zm00001eb0… cta1 chitinase … 4.40 4.21e- 3
## 6 Leaf Downregulated Zm00001eb1… <NA> Protein tr… -1.48 7.05e- 8
## 7 Leaf Downregulated Zm00001eb1… <NA> MFS domain… -1.03 2.04e- 6
## 8 Leaf Upregulated Zm00001eb1… nactf108 NAC-transc… 1.17 6.26e-10
## 9 Leaf Upregulated Zm00001eb3… <NA> Epimerase … 1.32 1.89e- 7
## 10 Leaf Upregulated Zm00001eb2… GRMZM2G1722… Clp R doma… 0.781 2.49e- 7
## 11 Leaf Upregulated Zm00001eb3… nye1 non-yellow… 0.765 5.24e- 7
## 12 Leaf Upregulated Zm00001eb2… <NA> ATP-depend… 0.934 1.08e- 5
## 13 Leaf Upregulated Zm00001eb2… nactf74 NAC-transc… 0.874 2.26e- 4
## 14 Leaf Upregulated Zm00001eb0… cta1 chitinase … 1.60 1.20e- 3
## 15 Leaf Upregulated Zm00001eb0… pco102923 Aspergillu… 1.18 1.75e- 3
SAG GWAS - Sekhon et al. 2019
Purpose: Load GWAS-identified senescence genes and
convert to v5 with manual corrections for problematic mappings.
sag_gwas_Sekhon2019 <- read.csv(
file.path(paths$data, "senescence_Sekhon2019.csv")
) %>%
left_join(v4_v5, by = c(gene_v4 = "v4"))
# Manual corrections for genes with mapping issues
to_correct <- c(
"Zm00001d045427" = "Zm00001eb377840",
"Zm00001d044747" = "Zm00001eb371660",
"Zm00001d039384" = "Zm00001eb120170"
)
sag_gwas_Sekhon2019$v5[
sag_gwas_Sekhon2019$gene_v4 %in% names(to_correct)
] <- to_correct
cat("Total SAG GWAS genes (Sekhon 2019):",
length(unique(sag_gwas_Sekhon2019$v5)), "\n")
## Total SAG GWAS genes (Sekhon 2019): 75
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% sag_gwas_Sekhon2019$v5
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 3 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Downregulated Zm00001eb41… <NA> <NA> -2.23 1.54e-3
## 2 -P Upregulated Zm00001eb28… gst41 glutathion… 3.02 3.54e-6
## 3 Leaf Downregulated Zm00001eb12… <NA> Beta-glucu… -0.938 6.43e-8
Stay-Green Network - Sekhon et al. 2019
Supplemental Data Set 5. Sekhon et al. 2019 https://academic.oup.com/plcell/article/31/9/1968/5985812
Purpose: Load genes in the stay-green regulatory
network and convert to v5.
staygreen_network <- read.csv(
file.path(paths$data, "staygreen_network_sekhon2019.csv"),
na.strings = ""
) %>%
janitor::clean_names() %>%
inner_join(v4_v5, by = c(gene = "v4"))
cat("Total stay-green network genes:",
length(unique(staygreen_network$v5)), "\n")
## Total stay-green network genes: 695
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% staygreen_network$v5
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 82 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Downregulated Zm00001eb06… <NA> GDSL ester… -2.83 2.50e-4
## 2 -P Downregulated Zm00001eb25… <NA> <NA> -2.02 1.07e-2
## 3 -P Downregulated Zm00001eb12… <NA> Carboxypep… -2.91 4.17e-2
## 4 -P Downregulated Zm00001eb13… phys2 phytase2 -2.20 4.97e-2
## 5 -P Upregulated Zm00001eb28… <NA> <NA> 3.39 1.58e-7
## 6 -P Upregulated Zm00001eb35… gst15 glutathion… 2.33 4.14e-6
## 7 -P Upregulated Zm00001eb30… chn16 chitinase16 2.88 4.29e-6
## 8 -P Upregulated Zm00001eb40… IDP132 <NA> 4.29 8.66e-6
## 9 -P Upregulated Zm00001eb06… IDP7286 Bowman-Bir… 3.55 9.55e-6
## 10 -P Upregulated Zm00001eb33… lrk1 Ser/Thr re… 2.08 3.83e-5
## # ℹ 72 more rows
Natural Senescence Genes - Zhang et al. 2014
Purpose: Load genes associated with natural
senescence processes and convert to v5.
natural_senescence <- read.csv(
file.path(paths$data, "natural_senescence.csv")
) %>%
janitor::clean_names() %>%
inner_join(v3_v5, by = c(gene_id = "v3"))
cat("Total natural senescence genes:",
length(unique(natural_senescence$v5)), "\n")
## Total natural senescence genes: 4722
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% natural_senescence$v5
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 629 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Downregulated Zm00001eb29… peamt2 phosphoeth… -6.80 0.0000172
## 2 -P Downregulated Zm00001eb16… rca3 RUBISCO ac… -3.69 0.0000740
## 3 -P Downregulated Zm00001eb20… <NA> Chlorophyl… -2.55 0.0000932
## 4 -P Downregulated Zm00001eb31… nactf14 NAC-transc… -2.57 0.000152
## 5 -P Downregulated Zm00001eb12… <NA> <NA> -3.17 0.000309
## 6 -P Downregulated Zm00001eb05… GRMZM2G1336… Integral m… -2.16 0.000648
## 7 -P Downregulated Zm00001eb04… <NA> Long-chain… -2.38 0.000757
## 8 -P Downregulated Zm00001eb01… tps8 terpene sy… -3.96 0.000999
## 9 -P Downregulated Zm00001eb03… <NA> MATH domai… -2.16 0.00141
## 10 -P Downregulated Zm00001eb23… <NA> <NA> -2.78 0.00177
## # ℹ 619 more rows
SAG Hub Genes - Ojeda et al. 2026
Purpose: Load hub genes from senescence
co-expression network.
sag_hub_ojeda2026 <- read.table(
file.path(paths$data, "senesence_hubs_ojeda2026.tab"),
header = TRUE,
sep = "\t",
fill = TRUE
)
cat("Total SAG hub genes (Ojeda 2026):",
length(unique(sag_hub_ojeda2026$gene)), "\n")
## Total SAG hub genes (Ojeda 2026): 63
sag_hubs <- sag_hub_ojeda2026 %>%
left_join(
v3_v5 %>%
separate_longer_delim(v5, delim = ",") %>%
dplyr::select(v3 = v3, v5) %>%
group_by(v5) %>%
dplyr::slice(1),
by = c(gene = "v5")
) %>%
left_join(
v4_v5 %>%
separate_longer_delim(v5, delim = ",") %>%
dplyr::select(v4 = v4, v5) %>%
group_by(v5) %>%
dplyr::slice(1),
by = c(gene = "v5")
) %>%
dplyr::select(v3,v4,gene,everything()) %>%
inner_join (effects_df %>%
filter(
is_hiconf_DEG,
gene %in% sag_hub_ojeda2026$gene
)
) %>%
select(predictor, regulation,v3,v4,gene,locus_label,symbol,functional_process, logFC, adj.P.Val,
desc_merged) %>%
ungroup() %>%
group_by(predictor, regulation) %>%
arrange(predictor, regulation, adj.P.Val)
write.csv(sag_hubs,
file.path(paths$intermediate, "sag_hub_ojeda2026_hiconf_degs.csv"),
row.names = FALSE)
nrow(sag_hubs)
## [1] 14
Photosynthesis Genes - Ojeda et al. 2026
Purpose: Load photosynthesis-related genes for
expression index calculation.
photosynthesis_ojeda2026 <- read.table(
file.path(paths$data, "photosynthesis_expression_index_genes_ojeda2026.txt"),
header = FALSE,
sep = "\t"
)
colnames(photosynthesis_ojeda2026) <- "gene"
cat("Total photosynthesis genes (Ojeda 2026):",
nrow(photosynthesis_ojeda2026), "\n")
## Total photosynthesis genes (Ojeda 2026): 172
photo_ojeda2026 <- effects_df %>%
filter(
is_hiconf_DEG,
gene %in% photosynthesis_ojeda2026$gene
) %>%
select(predictor, regulation, gene,locus_label, logFC, adj.P.Val,
desc_merged) %>%
ungroup() %>%
group_by(predictor, regulation) %>%
arrange(predictor, regulation, adj.P.Val)
write.csv(photo_ojeda2026,
file.path(paths$intermediate, "photo_ojeda2026_hiconf_degs.csv"),
row.names = FALSE)
nrow(photo_ojeda2026)
## [1] 6
GO-Based Gene Sets
Purpose: Extract photosynthesis and leaf senescence
genes from GO annotations.
GO terms used: - Photosynthesis: GO:0015979, GO:0019684, GO:0009768 - Leaf senescence: GO:0010150
photosynthesis_genes <- TERM2GENE %>%
filter(GO %in% c("GO:0015979", "GO:0019684", "GO:0009768")) %>%
pull(gene) %>%
unique()
leaf_senescence_genes <- TERM2GENE %>%
filter(GO %in% c("GO:0010150")) %>%
pull(gene) %>%
unique()
cat("GO photosynthesis genes:", length(photosynthesis_genes), "\n")
## GO photosynthesis genes: 672
cat("GO leaf senescence genes:", length(leaf_senescence_genes), "\n")
## GO leaf senescence genes: 200
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% photosynthesis_genes
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 21 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Downregulated Zm00001eb1… rca3 RUBISCO ac… -3.69 7.40e- 5
## 2 -P Downregulated Zm00001eb2… <NA> Chlorophyl… -2.55 9.32e- 5
## 3 -P Downregulated Zm00001eb3… lhcb10 light harv… -2.29 5.91e- 3
## 4 -P Upregulated Zm00001eb2… pfk1 phosphofru… 2.62 3.94e- 9
## 5 -P Upregulated Zm00001eb3… <NA> Phosphoeno… 5.36 6.50e- 9
## 6 -P Upregulated Zm00001eb2… <NA> ATP-depend… 2.23 2.33e- 4
## 7 -P Upregulated Zm00001eb2… glp1 germin-lik… 2.44 5.39e- 4
## 8 Inv4m Downregulated Zm00001eb3… Zm00001d008… Peptidyl-p… -2.80 2.63e- 6
## 9 Leaf Downregulated Zm00001eb0… dpr1 dihydrodip… -1.16 9.16e-12
## 10 Leaf Downregulated Zm00001eb0… pcr2 protochlor… -0.706 1.70e-11
## # ℹ 11 more rows
effects_df %>%
filter(
is_hiconf_DEG,
gene %in% leaf_senescence_genes
) %>%
select(predictor, regulation, gene:logFC, adj.P.Val,
desc_merged) %>%
group_by(predictor) %>%
arrange(predictor, regulation, adj.P.Val) %>%
tibble()
## # A tibble: 21 × 7
## predictor regulation gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 -P Downregulated Zm00001eb1… rca3 RUBISCO ac… -3.69 7.40e-5
## 2 -P Upregulated Zm00001eb1… <NA> <NA> 4.38 1.69e-7
## 3 -P Upregulated Zm00001eb1… <NA> Aspergillu… 7.50 9.10e-7
## 4 -P Upregulated Zm00001eb1… cl5103_-2b Methyleste… 2.25 4.07e-6
## 5 -P Upregulated Zm00001eb3… <NA> Methyleste… 2.65 1.09e-5
## 6 -P Upregulated Zm00001eb0… GRMZM2G1687… Aspergillu… 3.47 1.80e-4
## 7 -P Upregulated Zm00001eb2… nactf30 NAC-transc… 2.69 1.31e-3
## 8 -P Upregulated Zm00001eb0… dapat3 diaminopim… 3.73 2.54e-2
## 9 Leaf Downregulated Zm00001eb3… <NA> DIMBOA UDP… -0.737 9.52e-9
## 10 Leaf Downregulated Zm00001eb3… <NA> C3H1-type … -0.817 1.10e-7
## # ℹ 11 more rows
Expression Index Calculation
Purpose: Calculate normalized expression indices for
photosynthesis and senescence gene sets across samples.
Approach: Use mean expression of genes within each
set, normalize per set, and summarize by leaf stage and treatment.
# Define gene sets restricted to expression matrix universe
gene_set_list <- list(
photosynthesis = intersect(photosynthesis_ojeda2026$gene, universe),
senescence = intersect(sag_hub_ojeda2026$gene, universe)
)
cat("\n=== Gene Set Sizes ===\n")
##
## === Gene Set Sizes ===
print(sapply(gene_set_list, length))
## photosynthesis senescence
## 165 56
# Calculate raw expression indices
xpindex_data <- calculate_xpindex(
expression_matrix = normalized_expression,
gene_sets = gene_set_list,
method = "mean"
)
## photosynthesis: 165/165 genes found
## senescence: 56/56 genes found
# Normalize indices per set
xpindex_norm <- normalize_xpindex(
xpindex_data = xpindex_data,
index_names = names(gene_set_list),
method = "per_set"
)
# Summarize by leaf stage and treatment
xpindex_summary <- summarize_xpindex(
xpindex_data = xpindex_norm,
index_names = names(gene_set_list),
group_vars = c("leaf_stage", "treatment")
) %>%
mutate(
xpindex = factor(
tools::toTitleCase(xpindex),
levels = c("Photosynthesis", "Senescence"),
labels = c("Photosynthesis", "Leaf Senescence")
)
)
Pigment Biosynthesis Indices (CornCyc Pathways)
Purpose: Calculate xpindex for pigment biosynthesis
using CornCyc pathway annotations.
Approach: Extract gene sets from CornCyc, calculate
mean xpindex, normalize per_set, and visualize developmental
trajectories.
corncyc <- read.table(
file.path(paths$data, "corncyc_pathways.txt"),
sep = "\t",
header = TRUE,
quote = ""
)
pathway_ids <- list(
chlorophyll = c("CHLOROPHYLL-SYN", "PWY-5068","PWY-5188"),
carotenoid = "CAROTENOID-PWY",
anthocyanin = "PWY-5125",
flavonoids = "PWY1F-FLAVSYN"
)
corncyc_gene_sets <- lapply(names(pathway_ids), function(pigment) {
corncyc_genes <- corncyc %>%
filter(Pathway.id %in% pathway_ids[[pigment]]) %>%
pull(Protein.id) %>%
gsub("ZM00001EB", "Zm00001eb", .) %>%
gsub("_P\\d+$", "", ., perl = TRUE)
missing <- corncyc_genes[!grepl("Zm00001eb", corncyc_genes)]
if (length(missing) > 0) {
cat(pigment, "- Protein.id missing:",
paste(missing, collapse = ", "), "\n")
}
corncyc_genes[grepl("Zm00001eb", corncyc_genes)] %>% unique()
})
## chlorophyll - Protein.id missing: GDQC-106314-MONOMER, GDQC-112611-MONOMER, MONOMERDQC-5752, GDQC-114541-MONOMER, GDQC-108880-MONOMER
## carotenoid - Protein.id missing: MONOMER-15665, GBWI-61910-MONOMER, GBWI-61910-MONOMER, MONOMER-15665, MONOMER-15661
names(corncyc_gene_sets) <- names(pathway_ids)
{
cat("\n=== Pigment Gene Set Sizes ===\n")
print(sapply(corncyc_gene_sets, length))
}
##
## === Pigment Gene Set Sizes ===
## chlorophyll carotenoid anthocyanin flavonoids
## 40 35 15 51
corncyc_xpindex <- calculate_xpindex(
expression_matrix = normalized_expression,
gene_sets = corncyc_gene_sets,
method = "mean"
)
## chlorophyll: 31/40 genes found
## carotenoid: 29/35 genes found
## anthocyanin: 10/15 genes found
## flavonoids: 32/51 genes found
corncyc_xpindex_norm <- normalize_xpindex(
xpindex_data = corncyc_xpindex,
index_names = names(pathway_ids),
expression_matrix = normalized_expression,
gene_sets = corncyc_gene_sets,
method = "per_set"
)
{
cat("\n=== Normalized Pigment Index Ranges ===\n")
corncyc_xpindex_norm %>%
summarise(across(
all_of(names(pathway_ids)),
list(min = min, max = max),
.names = "{.col}_{.fn}"
)) %>%
pivot_longer(
cols = everything(),
names_to = c("index", "stat"),
names_pattern = "(.+)_(min|max)"
) %>%
pivot_wider(
names_from = stat,
values_from = value
) %>%
print()
}
##
## === Normalized Pigment Index Ranges ===
## # A tibble: 4 × 3
## index min max
## <chr> <dbl> <dbl>
## 1 chlorophyll 0 1
## 2 carotenoid 0 1
## 3 anthocyanin 0 1
## 4 flavonoids 0 1
corncyc_summary <- summarize_xpindex(
xpindex_data = corncyc_xpindex_norm,
index_names = names(pathway_ids),
group_vars = "leaf_stage"
)
effects_df %>%
filter(is_hiconf_DEG,
predictor=="Leaf",
gene %in% corncyc_gene_sets$chlorophyll) %>%
select(gene:logFC,adj.P.Val,desc_merged) %>% tibble()
## # A tibble: 0 × 5
## # ℹ 5 variables: gene <chr>, locus_symbol <chr>, desc_merged <chr>,
## # logFC <dbl>, adj.P.Val <dbl>
effects_df %>%
filter(gene =="Zm00001eb191650") %>%
select(gene:logFC,adj.P.Val,desc_merged) %>% tibble()
## # A tibble: 4 × 5
## gene locus_symbol desc_merged logFC adj.P.Val
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Zm00001eb191650 phos2 phosphate transporter2 0.533 0.000638
## 2 Zm00001eb191650 phos2 phosphate transporter2 -0.989 0.0829
## 3 Zm00001eb191650 phos2 phosphate transporter2 -0.457 0.556
## 4 Zm00001eb191650 phos2 phosphate transporter2 0.693 0.607
Pigment Plot
base_size <- 30
p_pigments <- corncyc_summary %>%
filter(xpindex != "flavonoids") %>%
mutate(xpindex = tools::toTitleCase(xpindex)) %>%
ggplot(aes(x = leaf_stage, y = mean, color = xpindex)) +
geom_line(linewidth = 1.5) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.2,
linewidth = 0.8
) +
labs(
x = "Leaf",
y = "Normalized Index",
color = "Pigment"
) +
scale_color_manual(
values = c(
"Chlorophyll" = "darkgreen",
"Carotenoid" = "gold",
"Anthocyanin" = "purple4"
)
) +
theme_classic(base_size = base_size) +
theme(legend.position = c(0.9, 0.9))
print(p_pigments)

ggsave(
file.path(paths$figures, "corncyc_pigment_biosynthesis_indices.png"),
plot = p_pigments,
width = 12,
height = 8,
dpi = 300
)
p_by_treatment <- ggplot(
xpindex_summary,
aes(x = leaf_stage, y = mean, color = xpindex, linetype = treatment)
) +
geom_line(linewidth = 1.5) +
geom_point(size = 4) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.2,
linewidth = 1
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen",
"Leaf Senescence" = "orange")
) +
scale_linetype_manual(
values = c("+P" = "solid", "-P" = "dashed")
) +
labs(
x = "Leaf",
y = "Normalized Expression Index",
color = NULL,
linetype = "Treatment"
) +
theme_classic(base_size = base_size) +
theme(legend.position = "right")
print(p_by_treatment)

Main Figure: Left Panel Components
Purpose: Create publication-quality composite figure
showing normalized xpindex, individual gene trajectories, and marker
gene validation.
Panel A: Normalized Gene Set Indices
xpindex_summary_pooled <- summarize_xpindex(
xpindex_data = xpindex_norm,
index_names = names(gene_set_list),
group_vars = "leaf_stage"
) %>%
mutate(
xpindex = factor(
tools::toTitleCase(xpindex),
levels = c("Photosynthesis", "Senescence"),
labels = c("Photosynthesis", "Leaf Senescence")
)
)
# --- Compute correlation stats ---
photo_cor <- cor.test(xpindex_norm$leaf_stage, xpindex_norm$photosynthesis, method = "pearson")
sen_cor <- cor.test(xpindex_norm$leaf_stage, xpindex_norm$senescence, method = "pearson")
# Format to 2 significant figures
format_2sig <- function(x) {
x_sig <- signif(x, digits = 2)
format(x_sig, scientific = FALSE, trim = TRUE, nsmall = 0)
}
# Create plotmath label as character string
make_label <- function(r_val, p_val) {
r_str <- format_2sig(r_val)
if (p_val < 0.01) {
exp_val <- floor(log10(p_val))
coeff <- round(p_val / 10^exp_val, 2)
if (coeff >= 10) {
coeff <- coeff / 10
exp_val <- exp_val + 1
}
coeff_str <- format_2sig(coeff)
paste0("italic(r) == ", r_str, " * ',' ~ italic(p) == ", coeff_str, " %*% 10^", exp_val)
} else {
p_str <- format_2sig(p_val)
paste0("italic(r) == ", r_str, " * ',' ~ italic(p) == ", p_str)
}
}
# Build annotation data frame
annotation_df <- data.frame(
x = 1.0,
y = c(0.51, 0.35), # adjust y-spacing here
label = c(
make_label(photo_cor$estimate, photo_cor$p.value),
make_label(sen_cor$estimate, sen_cor$p.value)
),
stringsAsFactors = FALSE
)
# --- Plot ---
base_family <- "Helvetica"
up <- 58
p_normalized <- ggplot(
xpindex_summary_pooled,
aes(x = leaf_stage, y = mean, color = xpindex)
) +
geom_line(linewidth = 1.5) +
geom_point(size = 4) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.2,
linewidth = 1
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen",
"Leaf Senescence" = "orange")
) +
# Single annotation layer using data frame
geom_text(
data = annotation_df,
aes(x = x, y = y, label = label),
hjust = 0,
vjust = 0,
size = 5,
parse = TRUE,
family = base_family,
inherit.aes = FALSE
) +
labs(
x = "Leaf",
y = "Normalized Expression Index",
color = NULL
) +
ggpubr::theme_classic2(base_size = base_size) +
theme(
legend.position = c(0.7, 0.9),
legend.text = element_text(size = 20),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(r = 50, unit = "pt"),
plot.margin = margin(5.5, 5.5, up, 7, "pt")
)
print(p_normalized)

Panel B: Spaghetti Plot (Individual Gene Trajectories)
xp_photo_genes <- intersect(
gene_set_list$photosynthesis,
rownames(normalized_expression)
)
xp_photosynthesis <- normalized_expression[xp_photo_genes, ] %>%
as.data.frame() %>%
tibble::rownames_to_column("gene") %>%
pivot_longer(-gene, names_to = "sample", values_to = "expression") %>%
mutate(
leaf_stage = as.numeric(sub(".*L([0-9]+)[LR]$", "\\1", sample))
) %>%
group_by(gene) %>%
mutate(centered_expr = expression - mean(expression)) %>%
ungroup() %>%
group_by(gene, leaf_stage) %>%
summarise(mean_expr = mean(centered_expr), .groups = "drop") %>%
arrange(gene)
xp_senescence_genes <- intersect(
gene_set_list$senescence,
rownames(normalized_expression)
)
xp_senescence <- normalized_expression[xp_senescence_genes, ] %>%
as.data.frame() %>%
tibble::rownames_to_column("gene") %>%
pivot_longer(-gene, names_to = "sample", values_to = "expression") %>%
mutate(
leaf_stage = as.numeric(sub(".*L([0-9]+)[LR]$", "\\1", sample))
) %>%
group_by(gene) %>%
mutate(centered_expr = expression - mean(expression)) %>%
ungroup() %>%
group_by(gene, leaf_stage) %>%
summarise(mean_expr = mean(centered_expr), .groups = "drop") %>%
arrange(gene)
{
cat("\n=== Spaghetti Plot Gene Counts ===\n")
cat("Photosynthesis genes:", length(xp_photo_genes), "\n")
cat("Senescence genes:", length(xp_senescence_genes), "\n")
}
##
## === Spaghetti Plot Gene Counts ===
## Photosynthesis genes: 165
## Senescence genes: 56
p_spaghetti <- xp_photosynthesis %>%
ggplot(aes(x = leaf_stage, y = mean_expr, group = gene)) +
geom_line(alpha = 0.1, linewidth = 0.5, color = "darkgreen") +
geom_smooth(
data = xp_photosynthesis,
aes(x = leaf_stage, y = mean_expr, group = 1),
method = "lm",
se = FALSE,
linewidth = 2,
color = "darkgreen"
) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 2) +
geom_line(
data = xp_senescence,
aes(x = leaf_stage, y = mean_expr, group = gene),
alpha = 0.1,
linewidth = 0.5,
color = "orange"
) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 2) +
geom_smooth(
data = xp_senescence,
aes(x = leaf_stage, y = mean_expr, group = 1),
method = "lm",
se = FALSE,
linewidth = 2,
color = "orange"
) %>%
with_shadow(colour = "black", x_offset = 0, y_offset = 0, sigma = 2) +
labs(
x = "Leaf",
y = expression("Centered " * log[2] * "(CPM)")
) +
coord_cartesian(ylim = c(-1, 1)) +
scale_y_continuous(position = "right") +
ggpubr::theme_classic2(base_size = base_size) +
theme(plot.margin = margin(5.5, 7, up, 5.5, "pt"))
print(p_spaghetti)

Combine Top Panels
p_top_panel <- ggarrange(
p_normalized + rremove("xlab"),
p_spaghetti + rremove("xlab"),
ncol = 2
) %>%
annotate_figure(
bottom = text_grob("Leaf", size = base_size, vjust = -2)
)
print(p_top_panel)

Panel C: Marker Gene Ilustration
marker_genes <- data.frame(
locus_label = c("pep1", "salt1",
"ssu1", "mir3",
"nye2", "sgrl1",
"cab1", "dapat3"),
gene_set = c("Photosynthesis", "Senescence",
"Photosynthesis", "Senescence",
"Senescence", "Senescence",
"Photosynthesis", "Senescence"),
pair = c("Pair1", "Pair1",
"Pair2", "Pair2",
"Pair3", "Pair3",
"Pair4", "Pair4"),
stringsAsFactors = FALSE
)
marker_genes <- data.frame(
locus_label = c("pep1", "salt1",
"ssu1", "mir3",
"nye2", "sgrl1",
"cab1", "nactf108"),
gene_set = c("Photosynthesis", "Senescence",
"Photosynthesis", "Senescence",
"Senescence", "Senescence",
"Photosynthesis", "Senescence"),
pair = c("Pair1", "Pair1",
"Pair2", "Pair2",
"Pair3", "Pair3",
"Pair4", "Pair4"),
stringsAsFactors = FALSE
)
gene_ids <- effects_df %>%
filter(locus_label %in% marker_genes$locus_label, !is.na(locus_label)) %>%
select(gene, locus_label, desc_merged) %>%
distinct() %>%
group_by(locus_label) %>%
slice(1) %>%
ungroup()
marker_data <- lapply(seq_len(nrow(gene_ids)), function(i) {
g <- gene_ids$gene[i]
data.frame(
sample = colnames(normalized_expression),
expression = normalized_expression[g, ],
gene = g
)
}) %>%
bind_rows() %>%
left_join(gene_ids, by = "gene") %>%
left_join(marker_genes, by = "locus_label") %>%
mutate(
leaf_stage = as.numeric(sub(".*L([0-9]+)[LR]$", "\\1", sample))
)
marker_summary <- marker_data %>%
group_by(gene, locus_label, gene_set, desc_merged, pair, leaf_stage) %>%
summarise(
mean_expr = mean(expression),
se_expr = sd(expression) / sqrt(n()),
.groups = "drop"
)
pair_max_expr <- marker_summary %>%
group_by(pair) %>%
summarise(max_expr = max(mean_expr), .groups = "drop") %>%
arrange(desc(max_expr))
marker_summary_colored <- marker_summary %>%
select(pair, locus_label, desc_merged, gene_set) %>%
distinct() %>%
mutate(
wrapped_text = stringr::str_wrap(desc_merged, width = 25),
with_br = gsub("\n", "<br>", wrapped_text),
with_color = case_when(
gene_set == "Photosynthesis" ~
paste0("<span style='color:darkgreen'>", with_br, "</span>"),
gene_set == "Senescence" ~
paste0("<span style='color:orange'>", with_br, "</span>"),
TRUE ~ with_br
),
sort_order = ifelse(gene_set == "Photosynthesis", 1, 2)
)
pair_labels <- marker_summary_colored %>%
arrange(pair, sort_order) %>%
group_by(pair) %>%
summarise(
facet_label = paste(with_color, collapse = "<br><br>"),
.groups = "drop"
) %>%
left_join(pair_max_expr, by = "pair") %>%
arrange(desc(max_expr))
marker_summary <- marker_summary %>%
left_join(
pair_labels %>% select(pair, facet_label, max_expr),
by = "pair"
) %>%
mutate(
facet_label = factor(facet_label, levels = pair_labels$facet_label)
)
p_markers <- ggplot(
marker_summary,
aes(x = leaf_stage, y = mean_expr, color = gene_set, group = gene)
) +
geom_line(linewidth = 1.5) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = mean_expr - se_expr, ymax = mean_expr + se_expr),
width = 0.2,
linewidth = 0.8
) +
geom_text(
data = marker_summary %>% filter(leaf_stage == 1),
aes(x = leaf_stage,
y = mean_expr,
label = locus_label,
color = gene_set),
nudge_x = -0.5,
hjust = 1,
fontface = "bold",
size = 5
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen", "Senescence" = "orange")
) +
scale_x_continuous(expand = expansion(add = c(1.5, 0.2)), breaks = 1:4) +
facet_wrap(
~ facet_label,
nrow = 2,
scales = "free_y"
) +
labs(
x = "Leaf",
y = expression(log[2] * "(CPM)"),
color = NULL
) +
theme_classic(base_size = base_size) +
theme(
strip.background = element_blank(),
strip.text = element_markdown(size = 15, hjust = 0, face = "bold"),
plot.margin = margin(-50, 80, 5.5, 5.5, "pt"),
axis.title.y = element_text(margin = margin(r = -10)),
panel.spacing.x = unit(0, "cm"),
panel.spacing.y = unit(0, "cm"),
legend.position = "none"
)
print(p_markers)

Combine All Panels: Left Panel
left_panel <- cowplot::plot_grid(
p_top_panel, p_markers,
ncol = 1,
align = 'h',
axis = 'lr'
)
ggsave(
file.path(paths$figures, "left_panel.png"),
plot = left_panel,
width = 9,
height = 15,
dpi = 150
)
print(left_panel)

Build custom chlorophyll gene sets —
Manually curate sets
# Load CornCyc again if not in environment (safe re-load)
corncyc <- read.table(
file.path(paths$data, "corncyc_pathways.txt"),
sep = "\t",
header = TRUE,
quote = ""
)
# Extract synthesis genes: CHLOROPHYLL-SYN + PWY-5068 + PWY-5188
chloro_synthesis_genes <- corncyc %>%
filter(Pathway.id %in% c("CHLOROPHYLL-SYN", "PWY-5068", "PWY-5188")) %>%
pull(Protein.id) %>%
gsub("ZM00001EB", "Zm00001eb", .) %>%
gsub("_P\\d+$", "", ., perl = TRUE) %>%
grep("Zm00001eb", ., value = TRUE) %>%
unique()
# Extract degradation genes from PWY-5098
chloro_degradation_base <- corncyc %>%
filter(Pathway.id == "PWY-5098") %>%
pull(Protein.id) %>%
gsub("ZM00001EB", "Zm00001eb", .) %>%
gsub("_P\\d+$", "", ., perl = TRUE) %>%
grep("Zm00001eb", ., value = TRUE) %>%
unique()
# Manually add your curated degradation genes (using maizegdb_id = rownames in expression matrix)
# From Wei et al 2025
manual_degradation_genes <- c(
"Zm00001eb307720", # ZmCHL
"Zm00001eb047780", # ZmCHL2
"Zm00001eb231810", # ZmPPH
"Zm00001eb103480", # ZmSGR
"Zm00001eb076680", # ZmSGRL
"Zm00001eb140370" # ZmPAO
)
chloro_degradation_genes <- unique(c(chloro_degradation_base, manual_degradation_genes))
# Create named list
chloro_gene_sets <- list(
chlorophyll_synthesis = chloro_synthesis_genes,
chlorophyll_degradation = chloro_degradation_genes
)
effects_df %>%
filter(#is_hiconf_DEG,
predictor=="Leaf:-P",
gene %in% chloro_degradation_genes) %>%
select(gene:logFC,adj.P.Val,desc_merged) %>% tibble()
## # A tibble: 0 × 5
## # ℹ 5 variables: gene <chr>, locus_symbol <chr>, desc_merged <chr>,
## # logFC <dbl>, adj.P.Val <dbl>
Calculate and normalize xpindex
cat("=== Chlorophyll Gene Sets ===\n")
## === Chlorophyll Gene Sets ===
cat("Synthesis genes:", length(chloro_gene_sets$chlorophyll_synthesis), "\n")
## Synthesis genes: 40
cat("Degradation genes:", length(chloro_gene_sets$chlorophyll_degradation), "\n")
## Degradation genes: 10
chloro_xpindex <- calculate_xpindex(
expression_matrix = normalized_expression,
gene_sets = chloro_gene_sets,
method = "mean"
)
## chlorophyll_synthesis: 31/40 genes found
## chlorophyll_degradation: 9/10 genes found
chloro_xpindex_norm <- normalize_xpindex(
xpindex_data = chloro_xpindex,
index_names = names(chloro_gene_sets),
method = "per_set"
)
lm(data=chloro_xpindex,
chlorophyll_degradation~leaf_stage*treatment ) %>%
summary()
##
## Call:
## lm(formula = chlorophyll_degradation ~ leaf_stage * treatment,
## data = chloro_xpindex)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40137 -0.09293 0.02659 0.05372 0.80915
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.38762 0.09641 55.884 < 2e-16 ***
## leaf_stage 0.06527 0.03520 1.854 0.07130 .
## treatment-P -0.17022 0.14450 -1.178 0.24595
## leaf_stage:treatment-P 0.15017 0.05368 2.797 0.00796 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1928 on 39 degrees of freedom
## Multiple R-squared: 0.5173, Adjusted R-squared: 0.4802
## F-statistic: 13.93 on 3 and 39 DF, p-value: 2.535e-06
lm(data=chloro_xpindex,
chlorophyll_synthesis~leaf_stage*treatment ) %>%
summary()
##
## Call:
## lm(formula = chlorophyll_synthesis ~ leaf_stage * treatment,
## data = chloro_xpindex)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32398 -0.07748 0.02540 0.09880 0.41586
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.36190 0.08724 61.464 < 2e-16 ***
## leaf_stage -0.15342 0.03185 -4.816 2.23e-05 ***
## treatment-P 0.24275 0.13075 1.857 0.0709 .
## leaf_stage:treatment-P -0.11258 0.04858 -2.318 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1745 on 39 degrees of freedom
## Multiple R-squared: 0.6606, Adjusted R-squared: 0.6345
## F-statistic: 25.3 on 3 and 39 DF, p-value: 2.949e-09
Summarize for plotting
# Pooled (all leaves, ignore treatment)
summary_pooled <- summarize_xpindex(
xpindex_data = chloro_xpindex_norm,
index_names = names(chloro_gene_sets),
group_vars = "leaf_stage"
) %>%
mutate(
xpindex = case_when(
xpindex == "chlorophyll_synthesis" ~ "Synthesis",
xpindex == "chlorophyll_degradation" ~ "Degradation"
),
xpindex = factor(xpindex, levels = c("Synthesis", "Degradation"))
)
# By treatment
summary_by_treatment <- summarize_xpindex(
xpindex_data = chloro_xpindex_norm,
index_names = names(chloro_gene_sets),
group_vars = c("leaf_stage", "treatment")
) %>%
mutate(
xpindex = case_when(
xpindex == "chlorophyll_synthesis" ~ "Synthesis",
xpindex == "chlorophyll_degradation" ~ "Degradation"
),
xpindex = factor(xpindex, levels = c("Synthesis", "Degradation"))
)
Plotting 1
# --- Correlation stats for p_normalized (Photosynthesis / Senescence) ---
photo_cor <- cor.test(
xpindex_norm$leaf_stage,
xpindex_norm$photosynthesis,
method = "pearson")
sen_cor <- cor.test(
xpindex_norm$leaf_stage,
xpindex_norm$senescence,
method = "pearson")
annotation_df1 <- data.frame(
x = 1.0,
y = c(0.45, 0.38),
label = c(
make_label(photo_cor$estimate, photo_cor$p.value),
make_label(sen_cor$estimate, sen_cor$p.value)
),
stringsAsFactors = FALSE
)
synth_cor <- cor.test(
chloro_xpindex_norm$leaf_stage,
chloro_xpindex_norm$chlorophyll_synthesis,
method = "pearson")
degrad_cor <- cor.test(
chloro_xpindex_norm$leaf_stage,
chloro_xpindex_norm$chlorophyll_degradation,
method = "pearson")
annotation_df2 <- data.frame(
x = 1.0,
y = c(0.45, 0.38),
label = c(
make_label(synth_cor$estimate, synth_cor$p.value),
make_label(degrad_cor$estimate, degrad_cor$p.value)
),
stringsAsFactors = FALSE
)
p_normalized <- ggplot(
xpindex_summary_pooled,
aes(x = leaf_stage, y = mean, color = xpindex)
) +
geom_line(linewidth = 1.5) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.2,
linewidth = 1
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen",
"Leaf Senescence" = "orange")
) +
geom_text(
data = annotation_df1,
aes(x = x, y = y, label = label),
hjust = 0,
vjust = 0,
size = 4,
parse = TRUE,
family = base_family,
inherit.aes = FALSE
) +
labs(
x = "Leaf",
y = "Normalized Expression",
color = NULL
) +
ggpubr::theme_classic2(base_size = base_size) +
theme(
legend.position = c(0.7, 0.9),
legend.text = element_text(size = 25),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(r = 50, unit = "pt"),
plot.margin = margin(5.5, 5.5, up, 7, "pt")
)
print(p_normalized)

# --- Plot 2: p_chloro_pooled ---
base_size <- 30
p_chloro_pooled <- summary_pooled %>%
ggplot(aes(x = leaf_stage, y = mean, color = xpindex)) +
geom_line(linewidth = 1.5) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
width = 0.2,
linewidth = 0.8
) +
scale_color_manual(
values = c("Synthesis" = "darkgreen",
"Degradation" = "darkorange")
) +
geom_text(
data = annotation_df2,
aes(x = x, y = y, label = label),
hjust = 0,
vjust = 0,
size = 4,
parse = TRUE,
family = base_family,
inherit.aes = FALSE
) +
labs(
y = "Normalized Expression",
x = "Leaf",
color = NULL
) +
theme_classic(base_size = base_size) +
theme(legend.position = "none")
print(p_chloro_pooled)

ggsave(
file.path(paths$figures, "chlorophyll_synthesis_vs_degradation_pooled.png"),
plot = p_chloro_pooled,
width = 10,
height = 6,
dpi = 300
)
Plotting 2 joins
pd <- position_dodge(width = 0.2)
# Plot 2: By treatment
p_chloro_treatment <- summary_by_treatment %>%
ggplot(aes(x = leaf_stage, y = mean,
color = xpindex,
fill = xpindex,
linetype = treatment,
shape=treatment)) +
geom_line(linewidth = 2) +
geom_point(position = pd,
size = 5) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
position = pd,
width = 0.2,
linewidth = 1,
linetype = "solid",
) +
scale_x_continuous(
expand=expansion(add = c(0.2,0.2))
) +
scale_y_continuous(
expand=expansion(add = c(0,0.2)),
breaks = 2*(0:5)/10
) +
scale_color_manual(
values = c("Synthesis" = "darkgreen",
"Degradation" = "darkorange"),
labels = c("Chlorophyll\nSynthesis","Degradation")
) +
scale_fill_manual(
values = c("Synthesis" = "darkgreen",
"Degradation" = "darkorange", guide= "none")
) +
scale_linetype_manual(
values = c("+P" = "solid", "-P" = "dashed", guide= "none")
) +
scale_shape_manual(
values = c("+P" = 19, "-P" = 25)
) +
labs(
x = "Leaf",
y = "Normalized Expression",
color = NULL,
shape = NULL,
) +
guides(
fill = "none",
linetype = "none",
shape = guide_legend(
override.aes = list(fill = "black"),
order = 2,
#label.position = "left",
#label.hjust = 1 ,
reverse = TRUE
),
color = guide_legend(
override.aes = list(shape = NA),
order = 1,
# label.position = "left",
# label.hjust = 1,
)
)+
theme_classic(base_size = base_size) +
theme(
legend.title = element_blank(),
legend.position = c(0.1, 1),
legend.justification = c(0, 1),
legend.text = element_text(
size = 25,
margin = margin(l =unit(2, "lines"))
),
legend.box = "horizontal",
legend.box.just = "bottom",
legend.margin = margin(0, 0, 0, 0),
legend.spacing.x = unit(1, "lines"),
legend.key.spacing.x = unit(2, "lines"),
legend.key.size = unit(2, "lines"),
legend.key = element_rect(fill = "transparent"),
legend.background = element_rect(fill = "transparent")
)
p_chloro <- ggpubr::ggarrange(p_chloro_pooled,
p_chloro_treatment,
align="hv")
p_chloro

ggsave(
file.path(paths$figures, "chlorophyll_synthesis_vs_degradation_by_treatment.png"),
plot = p_chloro,
width = 12,
height = 7,
dpi = 300
)
# --- Step 1: Prepare both datasets with consistent structure ---
p_photo_sene <- xpindex_summary %>%
ggplot(aes(x = leaf_stage, y = mean,
color = xpindex,
fill = xpindex,
linetype = treatment,
shape=treatment)) +
geom_line(linewidth = 1.5) +
geom_point(position = pd,
size = 5) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
position = pd,
width = 0.2,
linewidth = 1,
linetype = "solid",
) +
scale_y_continuous(
expand=expansion(add = c(0,0.2))
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen",
"Leaf Senescence" = "darkorange"),
) +
scale_fill_manual(
values = c("Photosynthesis" = "darkgreen",
"Leaf Senescence" = "darkorange"),
guide= "none"
) +
scale_linetype_manual(
values = c("+P" = "solid", "-P" = "dashed", guide= "none")
) +
scale_shape_manual(
values = c("+P" = 19, "-P" = 25),
guide="none"
) +
labs(
x = "Leaf",
y = "Normalized Expression",
color = NULL,
shape = NULL,
) +
guides(
fill = "none",
linetype = "none",
shape="none",
color = guide_legend(
override.aes = list(shape = NA),
order = 1,
#label.position = "left",
#label.hjust = 1,
)
)+
theme_classic(base_size = base_size) +
theme(
legend.title = element_blank(),
legend.position = c(0, 1),
legend.justification = c(0, 1),
legend.text = element_text(
size = 25,
margin = margin(l =unit(2, "lines")) ),
# legend.text.align = 1,
# legend.justification = "center",
# legend.box.just = "right",
# legend.margin = margin(0, 0, 0, 0),
# legend.spacing.y = unit(0, "lines"),
# legend.key.spacing.y = unit(0, "lines"),
legend.key.size = unit(2, "lines"),
legend.key = element_rect(fill = "transparent"),
legend.background = element_rect(fill = "transparent"),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank()
)
print(p_photo_sene)

p_combined_by_treatment <- ggpubr::ggarrange(
p_chloro_treatment + rremove("xlab"),
p_photo_sene + rremove("xlab"),
nrow=1,
widths= c(1.4,1)) %>%
annotate_figure(
bottom = text_grob("Leaf", size = base_size,
vjust = -0.25,
hjust = 0)
)
p_combined_by_treatment

ggsave(
file.path(paths$figures, "combined_by_treatment_xpindex.png"),
plot = p_combined_by_treatment,
width = 9.25,
height = 9,
dpi = 150
)
Pooled indices
# ------------------------------------------------------------
# 1. Formatting helpers (must be defined)
# ------------------------------------------------------------
# If not already defined, uncomment these:
format_2sig <- function(x) {
x_sig <- signif(x, digits = 2)
format(x_sig, scientific = FALSE, trim = TRUE, nsmall = 0)
}
make_p_label <- function(p_val) {
if (p_val < 0.01) {
exp_val <- floor(log10(p_val))
coeff <- round(p_val / 10^exp_val, 2)
if (coeff >= 10) {
coeff <- coeff / 10
exp_val <- exp_val + 1
}
coeff_str <- format_2sig(coeff)
paste0("italic(p) == ", coeff_str, " %*% 10^", exp_val)
} else {
p_str <- format_2sig(p_val)
paste0("italic(p) == ", p_str)
}
}
# ------------------------------------------------------------
# 2. Annotation data frames (one call each)
# ------------------------------------------------------------
annotation_left <- data.frame(
x = rep(2.5, 4),
y = c(0.65, 0.6, 0.2, 0.15),
label = c(
paste0("italic(r) == ", format_2sig(synth_cor$estimate)),
make_p_label(synth_cor$p.value),
paste0("italic(r) == ", format_2sig(degrad_cor$estimate)),
make_p_label(degrad_cor$p.value)
),
stringsAsFactors = FALSE
)
annotation_right <- data.frame(
x = rep(2.5, 4),
y = c(0.70, 0.65, 0.2, 0.15),
label = c(
paste0("italic(r) == ", format_2sig(photo_cor$estimate)),
make_p_label(photo_cor$p.value),
paste0("italic(r) == ", format_2sig(sen_cor$estimate)),
make_p_label(sen_cor$p.value)
),
stringsAsFactors = FALSE
)
# ------------------------------------------------------------
# 3. Common scale
# ------------------------------------------------------------
common_y_scale <- scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.2),
expand = expansion(add = c(0, 0.2))
)
size <- 6
# ------------------------------------------------------------
# 4. Left panel: Chlorophyll
# ------------------------------------------------------------
p_chloro_pooled_legend <- summary_pooled %>%
ggplot(aes(x = leaf_stage, y = mean, color = xpindex)) +
geom_line(linewidth = 1.5) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2, linewidth = 0.8) +
scale_color_manual(
values = c("Synthesis" = "darkgreen", "Degradation" = "darkorange"),
labels = c("Chlorophyll\nSynthesis", "Degradation")
) +
common_y_scale +
geom_text(
data = annotation_left,
aes(x = x, y = y, label = label),
hjust = 0,
vjust = 0,
size = size,
parse = TRUE,
family = base_family,
inherit.aes = FALSE
) +
labs(x = "Leaf", y = "Normalized Expression") +
theme_classic(base_size = base_size) +
theme(
legend.position = c(0, 1),
legend.justification = c(0, 1),
legend.text = element_text(
size = 25,
margin = margin(l = unit(2, "lines"))
),
legend.title = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.background = element_rect(fill = "transparent")
)
# ------------------------------------------------------------
# 5. Right panel: Photosynthesis & Senescence
# ------------------------------------------------------------
p_normalized_legend <- ggplot(xpindex_summary_pooled,
aes(x = leaf_stage, y = mean, color = xpindex)) +
geom_line(linewidth = 1.5) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2, linewidth = 1) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen", "Leaf Senescence" = "orange"),
labels = c("Photosynthesis", "Leaf Senescence")
) +
common_y_scale +
geom_text(
data = annotation_right,
aes(x = x, y = y, label = label),
hjust = 0,
vjust = 0,
size = size,
parse = TRUE,
family = base_family,
inherit.aes = FALSE
) +
labs(x = "Leaf", y = "Normalized Expression") +
ggpubr::theme_classic2(base_size = base_size) +
theme(
legend.position = c(0, 1),
legend.justification = c(0, 1),
legend.text = element_text(
size = 25,
margin = margin(l = unit(2, "lines"))
),
legend.title = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.background = element_rect(fill = "transparent"),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank()
)
# ------------------------------------------------------------
# 6. Combine and save
# ------------------------------------------------------------
p_combined_pooled_separate_legends <- ggpubr::ggarrange(
p_chloro_pooled_legend + rremove("xlab"),
p_normalized_legend + rremove("xlab"),
nrow = 1,
widths = c(1.33, 1)
) %>%
annotate_figure(
bottom = text_grob("Leaf", size = base_size, vjust = -0.3, hjust = 0)
)
ggsave(
file.path(paths$figures, "combined_pooled_xpindex_separate_legends.png"),
plot = p_combined_pooled_separate_legends,
width = 8.2,
height = 7,
dpi = 150
)
print(p_combined_pooled_separate_legends)
## Panel C: Marker Gene Ilustration with black facet labels
marker_genes <- data.frame(
locus_label = c("pep1", "salt1",
"ssu1", "mir3",
"nye2", "sgrl1",
"cab1", "nactf108"),
gene_set = c("Photosynthesis", "Senescence",
"Photosynthesis", "Senescence",
"Senescence", "Senescence",
"Photosynthesis", "Senescence"),
pair = c("Pair1", "Pair1",
"Pair2", "Pair2",
"Pair3", "Pair3",
"Pair4", "Pair4"),
stringsAsFactors = FALSE
)
# Trim ", chloroplastic" from desc_merged in effects_df before joining
gene_ids <- effects_df %>%
mutate(desc_merged = stringr::str_remove(desc_merged, ", chloroplastic") %>% # <-- (1) TRIM
stringr::str_remove( "Protein ") %>% # <-- (1) TRIM
stringr::str_remove(" protein")
) %>%
filter(locus_label %in% marker_genes$locus_label, !is.na(locus_label)) %>%
select(gene, locus_label, desc_merged) %>%
distinct() %>%
group_by(locus_label) %>%
slice(1) %>%
ungroup()
marker_data <- lapply(seq_len(nrow(gene_ids)), function(i) {
g <- gene_ids$gene[i]
data.frame(
sample = colnames(normalized_expression),
expression = normalized_expression[g, ],
gene = g
)
}) %>%
bind_rows() %>%
left_join(gene_ids, by = "gene") %>%
left_join(marker_genes, by = "locus_label") %>%
mutate(
leaf_stage = as.numeric(sub(".*L([0-9]+)[LR]$", "\\1", sample))
)
marker_summary <- marker_data %>%
group_by(gene, locus_label, gene_set, desc_merged, pair, leaf_stage) %>%
summarise(
mean_expr = mean(expression),
se_expr = sd(expression) / sqrt(n()),
.groups = "drop"
)
# Build PLAIN TEXT facet labels with line breaks and wrapped descriptions
pair_labels <- marker_summary %>%
distinct(pair, locus_label, gene, desc_merged, gene_set) %>%
mutate(
facet_label = paste0(locus_label, "\n", gene),
# Wrap ONLY the description part (e.g., width = 30)
# wrapped_desc = stringr::str_wrap(desc_merged, width = 30),
# Construct label: "locus gene\nwrapped description"
# facet_label = paste0(locus_label, " ", gene, "\n", wrapped_desc),
# trunc_desc = stringr::str_trunc(desc_merged, width = 27),
# facet_label = paste0(locus_label, " ", gene, "\n", trunc_desc),
order_key = ifelse(gene_set == "Photosynthesis", 1, 2)
) %>%
arrange(pair, order_key) %>%
group_by(pair) %>%
summarise(
facet_label = paste(facet_label, collapse = "\n\n"), # double line break between gene entries
.groups = "drop"
)
# Get max expression per pair for ordering
pair_max_expr <- marker_summary %>%
group_by(pair) %>%
summarise(max_expr = max(mean_expr), .groups = "drop") %>%
arrange(desc(max_expr))
# Join and set factor levels in correct order
pair_labels <- pair_labels %>%
left_join(pair_max_expr, by = "pair") %>%
arrange(desc(max_expr))
marker_summary <- marker_summary %>%
left_join(pair_labels %>% select(pair, facet_label, max_expr), by = "pair") %>%
mutate(
facet_label = factor(facet_label, levels = pair_labels$facet_label)
)
p_markers <- ggplot(
marker_summary,
aes(x = leaf_stage, y = mean_expr, color = gene_set, group = gene)
) +
geom_line(linewidth = 1.5) +
geom_point(size = 5) +
geom_errorbar(
aes(ymin = mean_expr - se_expr, ymax = mean_expr + se_expr),
width = 0.2,
linewidth = 0.8
) +
geom_text(
data = marker_summary %>% filter(leaf_stage == 1),
aes(x = leaf_stage,
y = mean_expr,
label = locus_label,
color = gene_set),
nudge_x = -0.25,
hjust = 1,
fontface = "bold.italic",
size = 5
) +
scale_color_manual(
values = c("Photosynthesis" = "darkgreen", "Senescence" = "orange")
) +
scale_x_continuous(expand = expansion(add = c(1.5, 0.2)), breaks = 1:4) +
facet_wrap(
~ facet_label,
nrow = 2,
scales = "free_y"
) +
labs(
x = "Leaf",
y = expression(log[2] * "(CPM)"),
color = NULL
) +
theme_classic(base_size = base_size) +
theme(
strip.background = element_blank(),
# (3) Use plain element_text with black color
strip.text = element_text(size = 16, hjust = 0, face = "bold", color = "black"),
plot.margin = margin(-50, 80, 5.5, 5.5, "pt"),
axis.title.y = element_text(margin = margin(r = -10)),
panel.spacing.x = unit(0, "cm"),
panel.spacing.y = unit(0, "cm"),
legend.position = "none"
)
print(p_markers)
## Print left panel 2
left_panel <- cowplot::plot_grid(
p_combined_pooled_separate_legends +
theme(plot.margin = margin(5.5, 80, up, 5.5, "pt")),
p_markers,
ncol = 1,
align = 'h',
axis = 'lr'
)
ggsave(
file.path(paths$figures, "left_panel2.svg"),
fix_text_size = FALSE,
plot = left_panel,
width = 9,
height = 15,
dpi = 150
)
print(left_panel)
## Caluclate model stats
# --- Define combined gene sets ---
panel_sets <- c(chloro_gene_sets, gene_set_list)
# --- Calculate and normalize expression indices ---
panel_data <- calculate_xpindex(
expression_matrix = normalized_expression,
gene_sets = panel_sets,
method = "mean"
)
## chlorophyll_synthesis: 31/40 genes found
## chlorophyll_degradation: 9/10 genes found
## photosynthesis: 165/165 genes found
## senescence: 56/56 genes found
panel_norm <- normalize_xpindex(
xpindex_data = panel_data,
index_names = names(panel_sets),
method = "per_set"
)
# --- Model 1: Full model with interaction (leaf_stage * treatment) ---
cat("=== MODEL 1: Leaf Stage × Treatment Interaction ===\n")
## === MODEL 1: Leaf Stage × Treatment Interaction ===
panel_models_interaction <- fit_xpindex_models(
xpindex_data = panel_norm,
index_names = names(panel_sets),
formula_rhs = "leaf_stage * treatment"
)
stats_interaction <- extract_xpindex_stats(panel_models_interaction)
print(stats_interaction)
## # A tibble: 12 × 8
## xpindex predictor effect std_err p_value r2 p_adj p_signif
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Chlorophyll_degrada… leaf_sta… 0.0376 0.0203 7.13e-2 0.517 8.56e-2 ns
## 2 Chlorophyll_degrada… leaf_sta… 0.0864 0.0309 7.96e-3 0.517 1.37e-2 *
## 3 Chlorophyll_degrada… treatmen… -0.0980 0.0832 2.46e-1 0.517 2.46e-1 ns
## 4 Chlorophyll_synthes… leaf_sta… -0.122 0.0253 2.23e-5 0.661 1.34e-4 ***
## 5 Chlorophyll_synthes… leaf_sta… -0.0893 0.0385 2.58e-2 0.661 3.87e-2 *
## 6 Chlorophyll_synthes… treatmen… 0.193 0.104 7.09e-2 0.661 8.56e-2 ns
## 7 Photosynthesis leaf_sta… -0.0812 0.0203 2.67e-4 0.700 1.07e-3 **
## 8 Photosynthesis leaf_sta… -0.120 0.0309 3.81e-4 0.700 1.14e-3 **
## 9 Photosynthesis treatmen… 0.305 0.0831 7.15e-4 0.700 1.72e-3 **
## 10 Senescence leaf_sta… 0.0930 0.0189 1.60e-5 0.715 1.34e-4 ***
## 11 Senescence leaf_sta… 0.0867 0.0288 4.57e-3 0.715 9.15e-3 **
## 12 Senescence treatmen… -0.129 0.0775 1.03e-1 0.715 1.13e-1 ns
# --- Model 2: Additive model (leaf_stage + treatment) ---
cat("\n=== MODEL 2: Additive Effects (Leaf Stage + Treatment) ===\n")
##
## === MODEL 2: Additive Effects (Leaf Stage + Treatment) ===
panel_models_additive <- fit_xpindex_models(
xpindex_data = panel_norm,
index_names = names(panel_sets),
formula_rhs = "leaf_stage + treatment"
)
stats_additive <- extract_xpindex_stats(panel_models_additive)
print(stats_additive)
## # A tibble: 8 × 8
## xpindex predictor effect std_err p_value r2 p_adj p_signif
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Chlorophyll_degrada… leaf_sta… 0.0747 0.0166 5.47e- 5 0.420 1.09e-4 ***
## 2 Chlorophyll_degrada… treatmen… 0.114 0.0369 3.58e- 3 0.420 5.72e-3 **
## 3 Chlorophyll_synthes… leaf_sta… -0.160 0.0201 8.73e-10 0.614 3.49e-9 ****
## 4 Chlorophyll_synthes… treatmen… -0.0267 0.0448 5.55e- 1 0.614 6.34e-1 ns
## 5 Photosynthesis leaf_sta… -0.133 0.0178 4.22e- 9 0.584 1.12e-8 ****
## 6 Photosynthesis treatmen… 0.0106 0.0397 7.91e- 1 0.584 7.91e-1 ns
## 7 Senescence leaf_sta… 0.130 0.0156 2.81e-10 0.649 2.25e-9 ****
## 8 Senescence treatmen… 0.0835 0.0348 2.13e- 2 0.649 2.84e-2 *
# --- Model 3: Leaf stage only (pooled across treatments) ---
cat("\n=== MODEL 3: Leaf Stage Only (Pooled) ===\n")
##
## === MODEL 3: Leaf Stage Only (Pooled) ===
panel_models_pooled <- fit_xpindex_models(
xpindex_data = panel_norm,
index_names = names(panel_sets),
formula_rhs = "leaf_stage"
)
stats_pooled <- extract_xpindex_stats(panel_models_pooled)
print(stats_pooled)
## # A tibble: 4 × 8
## xpindex predictor effect std_err p_value r2 p_adj p_signif
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Chlorophyll_degrada… leaf_sta… 0.0729 0.0182 2.52e- 4 0.282 2.52e-4 ***
## 2 Chlorophyll_synthes… leaf_sta… -0.160 0.0199 6.32e-10 0.610 2.40e-9 ****
## 3 Photosynthesis leaf_sta… -0.133 0.0176 2.63e- 9 0.583 3.51e-9 ****
## 4 Senescence leaf_sta… 0.129 0.0165 1.20e- 9 0.598 2.40e-9 ****
# --- Optional: Compare models via AIC (for one example index) ---
# Example: Chlorophyll synthesis
example_index <- "chlorophyll_synthesis"
if (example_index %in% names(panel_sets)) {
m_int <- panel_models_interaction[[example_index]]
m_add <- panel_models_additive[[example_index]]
m_pooled <- panel_models_pooled[[example_index]]
cat("\n=== AIC Comparison (", example_index, ") ===\n", sep = "")
AIC_comparison <- AIC(m_int, m_add, m_pooled)
rownames(AIC_comparison) <- c("Interaction", "Additive", "Pooled")
print(AIC_comparison)
}
##
## === AIC Comparison (chlorophyll_synthesis) ===
## df AIC
## Interaction 5 -42.23264
## Additive 4 -38.68468
## Pooled 3 -40.30490
# --- Export results (optional) ---
# write.csv(stats_interaction, file.path(paths$intermediate, "xpindex_interaction_stats.csv"), row.names = FALSE)
# write.csv(stats_pooled, file.path(paths$intermediate, "xpindex_pooled_stats.csv"), row.names = FALSE)